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ABSTRACT 

We present a simple framework for the growth and evolution of supermassive black holes (SMBHs) in the 
hierarchical structure formation paradigm, adopting the general idea that quasar activity is triggered in major 
mergers. In our model, black hole accretion is triggered during major mergers (mass ratio > 0.3) between host 
dark matter halos. The successive evolution of quasar luminosities follows a universal light curve form, during 
which the growth of the SMBH is modeled self-consistently: an initial exponential growth at constant Edding- 
ton ratio of order unity until it reaches the peak luminosity, followed by a power-law decay. Assuming that 
the peak luminosity correlates with the post-merger halo mass, we convolve the light curve with the triggering 
rate of quasar activity to predict the quasar luminosity function (LF). Our model reproduces the observed LF 
at 0.5 < z < 4.5 for the full luminosity ranges probed by current optical and X-ray surveys. At z < 0.5, our 
model underestimates the LF at Lboi < 10"*^ ergs"', allowing room for AGN activity triggered by secular pro- 
cesses instead of major mergers. At z > 4.5, in order to reproduce the observed quasar abundance, the typical 
quasar hosts must shift to lower mass halos, and/or minor mergers can also trigger quasar activity. Our model 
reproduces both the observed redshift evolution and luminosity dependence of the linear bias of quasar/ AGN 
clustering. Due to the scatter between instantaneous luminosity and halo mass, quasar/ AGN clustering weakly 
depends on luminosity at low to intermediate luminosities; but the linear bias rises rapidly with luminosity at 
the high luminosity end and at high redshift. In our model, the Eddington ratio distribution is roughly log- 
normal, which broadens and shifts to lower mean values from high luminosity quasars (Lboi ^ 10"*^ ergs"') to 
low luminosity AGNs (Lboi ^ 10"*^ ergs"'), in good agreement with observations. The model predicts that the 
vast majority of > 10^'^ M0 SMBHs were already in place by z = 1, and < 50% of them were in place by z = 2; 
but the less massive (< 10^ Mq) SMBHs were assembled more recently, likely more through secular processes 
than by major mergers - in accordance with the downsizing picture of SMBH assembly since the peak of bright 
quasar activity around z 2 - 3. 

Subject headings: black hole physics - galaxies: active - cosmology: observations - large-scale structure of 
universe - quasars: general - surveys 



1. INTRODUCTION 

The cosmic assembly and evolution of supermassive black 
holes (SMBHs) is a central topic in modern cosmology: in the 
local universe, SMBHs reside in almost every bulge dominant 
galaxy (e.g., Kormendy & Richstone 1995; Richstone et al. 
1998); and they likely played important roles during their 
coevolution with galaxy bulges (e.g.. Silk & Rees 1998; 
Wyithe&Loeb 2002, 2003; King 2003; Di Matteo et al. 
2005; Croton et al. 2006), as inferred from observed scaling 
relations between bulge properties and the mass of the central 
BH (e.g., Magorrian et al. 1998; Ferrarese & Merritt 2000; 
Gebhardt et al. 2000; Graham etal. 2001; Tremaine et al. 
2002; Marconi & Hunt 2003). It is also generally ac- 
cepted that the local dormant SMBH population was 
largely assembled via gas accretion during a luminous QSO 
phase' (e.g., Salpeter 1964; Zel'dovich & Novikov 1964; 
Lynden-Bell 1969; Soltan 1982; Small & Blandford 1992; 
Saluccietal. 1999; Yu & Tremaine 2002; Shankaretal. 
2004; Marconi et al. 2004). The statistics of the local dor- 
mant SMBHs and distant active QSOs therefore provide clues 
to the build-up of the local BH density across cosmic time, 

' In this paper we use the term "QSO" to refer to all actively accreting 
SMBHs. We adopt the convention that L|,oi = 10'*' ergs"' is the dividing 
hne between quasars and AGNs. This division is slightly lower than the 
traditional Seyfert/quasar division (Schmidt & Green 1983) of Mb = —23 or 
Lb„i-10l2Lo. 



and have implications for the hierarchical structure formation 
paradigm, as well as for the interactions between accreting 
SMBHs and galaxy bulges. 

One of the leading hypotheses for the QSO triggering 
mechanism is galaxy mergers (e.g., Hernquist 1989; Carlberg 
1990; Kauffmann & Haehnelt 2000; Hopkins et al. 2008, and 
references therein). In the hierarchical CDM paradigm, small 
structures merge to form large structures, and the merger 
rate of dark matter halos peaks at early times, in broad con- 
sistency with the observed peak of bright quasar activities. 
Gas-rich major mergers (i.e., those with comparable mass 
ratios) between galaxies provide an efficient way to chan- 
nel a large amount of gas into the central region to trig- 
ger starbursts and possibly feed rapid black hole growth. 
Observationally this merger hypothesis is supported by the 
ULIRG-quasar connection (e.g., Sanders & Mirabel 1996; 
Canalizo & Stockton 2001), the signature of recent mergers 
in quasar hosts (e.g., Bennert et al. 2008), and the small- 
scale overdensities of galaxies around luminous quasars or 
quasar binaries (e.g.. Fisher etal. 1996; Bahcall et al. 1997; 
Serberetal. 2006; Hennawi et al. 2006; Myers et al. 2008; 
Strand et al. 2008). These observations do not necessarily 
prove that mergers are directly responsible for triggering QSO 
activity, but they at least suggest that QSO activity is coinci- 
dent with mergers in many cases. 

The last decade has seen a number of models for QSO evo- 
lution based on the merger hypothesis (e.g., Haiman & Loeb 
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1998; Kauffmann & Haehnelt 2000; Wyithe & Loeb 2002, 
2003; Volonteri et al. 2003; Hopkins et al. 2008), which agree 
with the bulk of QSO observations reasonably well. Moti- 
vated by the success of these merger-based QSO models, we 
revisit this problem in this paper with improved observational 
data on SMBHs and QSOs from dedicated large surveys, and 
with updated knowledge of the merger rate as inferred from 
recent numerical simulations. The present work is different 
from previous studies in many aspects in both methodology 
and the observational tests used. Our goal is to check if a 
simple merger-based cosmological QSO framework can re- 
produce all the observed statistics of SMBHs (both active and 
dormant), and to put constraints on the physical properties 
of SMBH growth. Our model is observationally motivated, 
therefore we do not restrict ourselves to theoretical predic- 
tions of SMBH/QSO properties from either cosmological or 
merger event simulations. 

The paper is organized as follows. In §1.1 we review some 
aspects of local SMBH demographics and QSO luminosity 
function; in § 1 .2 we briefly review the current status of quasar 
clustering observations; §§1.3 and 1.4 discuss the halo and 
subhalo merger rate from numerical simulations. Our model 
formulism is presented in §2 and we present our fiducial 
model and compare with observations in §3. We discuss ad- 
ditional aspects of our model in §4 and conclude in §5. 

We use friends-of-friends (fof) mass with a link length 
b = 0.2 to define the mass of a halo (roughly corresponding to 
spherical overdensity mass with A ~ 200 times the mean mat- 
ter density), since the fitting formulae for the halo mass func- 
tion are nearly universal with this mass definition (Sheth et al. 
2001; Jenkins etal. 2001; Warren et al. 2006; Tinker etal. 
2008). For simplicity we neglect the slight difference be- 
tween /o/ mass and virial mass (Appendix A; and see White 
2002, for discussions on different mass definitions). Through- 
out the paper, L always refers to the bolometric luminosity, 
and we use subscripts b or x to denote B-band or X-ray lumi- 
nosity when needed. We adopt a flat ACDM cosmology with 
Qq = 0.26, Qa = 0.74, nb = 0.044, /; = 0.7, erg = 0.78, «, = 0.95. 
We use the Eisenstein & Hu (1999) transfer function to com- 
pute the linear power spectrum, and use the fitting formulae 
for halo abundance from Sheth & Tormen (1999) and for halo 
linear bias from Sheth et al. (2001) based on the ellipsoidal 
collapse model and tested against numerical simulations. 

1.1. Supermassive Black Hole Demographics 

In the local universe, the dormant BH mass function 
(BHMF) can be estimated by combining scaling relations be- 
tween BH mass and galaxy bulge properties, such as the M, - 
a relation (e.g., Gebhardt et al. 2000; Ferrarese & Merritt 
2000; Tremaine et al. 2002) or the M, - Lsph(Msph) relation 
(e.g., Magorrian et al. 1998; Marconi & Hunt 2003), with 
bulge luminosity or stellar mass/velocity dispersion functions 
(e.g., Salucci et al. 1999; Yu & Tremaine 2002; Marconi et al. 
2004; Shankar et al. 2004; McLure & Dunlop 2004; Yu & Lu 
2004, 2008; Tundo et al. 2007; Shankar et al. 2009b). Al- 
though some uncertainties exist on the usage of these scal- 
ing relations at both the high and low BH mass end (e.g., 
Lauer et al. 2007; Tundo et al. 2007; Hu 2008; Greene et al. 
2008; Graham 2008; Graham & Li 2009), the total BH mass 
density is estimated to be p, « 4 x lO^MoMpc""^ with an un- 
certainty of a factor ~ 1.5 (e.g. Yu & Lu 2008; Shankar et al. 
2009b) - but the shape of the local BHMF is more uncertain 
(see the discussion in Shankar et al. 2009b). 

It is now generally accepted that local SMBHs were once 



luminous QSOs (Salpeter 1964; Zel'dovich & Novikov 1964; 
Lynden-Bell 1969). An elegant argument tying the active 
QSO population to the local dormant SMBH population was 
proposed by Soltan (1982): if SMBHs grow mainly through a 
luminous QSO phase, then the accreted luminosity density of 
QSOs to z = should equal the local BH mass density: 

[°° dt f°° (1 -e)L 

/9.,acc=/ ^ 2^$(L,z)fl'L«p. , (1) 

Jo Jo 

where <l>(L,z) is the bolometric luminosity function (LF) perL 
interval. Given the observed QSO luminosity function, a rea- 
sonably good match between p. acc and p, can be achieved if 
the average radiative efficiency e ~ 0.1 (e.g., Yu & Tremaine 
2002; Shankar et al. 2004; Marconi et al. 2004; Hopkins et al. 
2007; Shankar et al. 2009b). The exact value of e, however, 
is subject to some uncertainties from the luminosity function 
and local BH mass density determination. 

Assuming that SMBH growth is through gas accretion, an 
extended version of the Soltan argument can be derived using 
the continuity equation- (cf. Small & Blandford 1992): 

dn(M„t) ^ d[n(M„t){M.)] _^ 

dt dM, " ' 

where n{M,,t) is the BH mass function per dM,, and (M.) 
is the mean accretion rate of BHs with mass {M,,M,+dM,) 
averaged over both active and inactive BHs. Given the lu- 
minosity function, and a model connecting luminosity to BH 
mass, one can derive the BH mass function at all times by 
integrating the continuity equation (e.g., Marconi et al. 2004; 
Merloni 2004; Shankar et al. 2009b). 

Since the local BHMF and the QSO luminosity func- 
tion together determine the assembly history of the cosmic 
SMBH population, any cosmological model of AGN/quasars 
must first reproduce the observed luminosity function 
(e.g., Haiman & Loeb 1998; Kauffmann & Haehnelt 2000; 
Wyithe & Loeb 2002, 2003; Volonteri et al. 2003; Lapi et al. 
2006; Hopkins et al. 2008; Shankar et al. 2009b). This is also 
one of the central themes of this paper 

There has been great progress in the measurements of 
the QSO luminosity function over wide redshift and lu- 
minosity ranges for the past decade, mostly in the op- 
tical band (e.g.. Fan et al. 2001, 2004; Wolf etal. 2003; 
Croometal. 2004; Richards et al. 2005, 2006; Jiang et al. 
2006; Fontanot et al. 2007), and the soft/hard X-ray band 
(e.g., Ueda et al. 2003; Hasinger et al. 2005; Silverman et al. 
2005; Bargeretal. 2005). While wide-field optical sur- 
veys provide the best constraints on the bright end of the 
QSO luminosity functions, deep X-ray surveys can probe 
the obscured faint end AGN population. Although it has 
been known for a while that the spatial density of optically- 
selected bright quasars peaks at redshift z ~ 2-3, the spa- 
tial density of fainter AGNs selected in the X-ray seems to 
peak at lower redshift (e.g., Steffen et al. 2003; Ueda etal. 
2003; Hasinger et al. 2005), a trend now confirmed in other 
observational bands as well (e.g., Bongiorno et al. 2007; 
Hopkins et al. 2007, and references therein) - the manifesta- 
tion of the so-called downsizing of the cosmic SMBH assem- 
bly. 

- Note that here the source term, i.e., the creation of a BH with mass 
M, by mechanisms other than gas accretion, is generally neglected (e.g.. 
Small & Blandford 1992; Merloni 2004; Shankar et al. 2009b). One possi- 
ble such mechanism is BH coalescence, which modifies the shape of the BH 
mass function but does not change the total BH mass density in the classical 
case (i.e., neglecting mass loss from gravitational radiation). Dry mergers at 
low redshift may act as a surrogate of shaping the BH mass function. 
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Hopkins et al. (2007) compiled QSO luminosity function 
data from surveys in various bands (optical. X-ray, IR, 
etc). Assuming a general AGN spectral energy distribu- 
tion (SED) and a column density distribution for obscuration, 
Hopkins et al. (2007) were able to fit a universal bolometric 
luminosity function based on these data. We adopt their com- 
piled bolometric LP data in our modeling. The advantage of 
using the bolometric LP data is that both unobscured and ob- 
scured SMBH growth are counted in the model; but we still 
suffer from the systematics of bolometric corrections. The 
nominal statistical/systematic uncertainty level of the bolo- 
metric LP is - 20-30% (cf., Shankar et al. 2009b). 

1.2. Quasar Clustering 

A new ingredient in SMBH studies, due to dedicated large- 
scale optical surveys such as the 2QZ (Croom et al. 2004) and 
SDSS (York et al. 2000), is quasar clustering. While quasar 
clustering studies can be traced back to more than two decades 
ago (e.g.. Shaver 1984), statistically significant results only 
came very recently (e.g., Porciani et al. 2004; Croom etal. 
2005; Porciani & Norberg 2006; Myers et al. 2006, 2007a,b; 
Shen et al. 2007; da Angela et al. 2008; Padmanabhan et al. 
2008; Shen et al. 2008b, 2009; Ross et al. 2009). 

Within the biased halo clustering picture (e.g.. Kaiser 
1984; Bardeen et al. 1986; Mo & White 1996), the observed 
QSO two-point correlation function implies that quasars 
live in massive dark matter halos and are biased tracers 
of the underlying dark matter. Por optically luminous 
quasars (Lboi > lO'^^ergs"'), the host halo mass inferred 
from clustering analysis is a few times 1O'^/2~'M0. Thus 
quasar clustering provides independent constraints on how 
SMBHs occupy dark matter halos, where the abundance 
and evolution of the latter population can be well stud- 
ied in analytical theories and numerical simulations (e.g.. 
Press & Schechter 1974; Bond et al. 1991; Lacey & Cole 
1993; Mo & White 1996; Sheth & Tormen 1999; Sheth et al. 
2001; Springel et al. 2005b). By comparing the relative abun- 
dance of quasars and their host halos, one can constrain the 
average quasar duty cycles or lifetimes (e.g.. Cole & Kaiser 
1989; Martini & Weinberg 2001; Haiman & Hui 2001) to be 
< 10'^ yr, which means bright quasars are short lived. 

More useful constraints on quasar models come from stud- 
ies of quasar clustering as a function of redshift and lu- 
minosity. The redshift evolution of quasar clustering con- 
strains the evolution of duty cycles of active accretion, while 
the luminosity dependence of quasar clustering puts con- 
straints on quasar light curves (LC). Successful cosmolog- 
ical quasar models therefore not only need to account for 
quasar abundances (LP), but also need to explain quasar clus- 
tering properties, both as function of redshift and luminos- 
ity. We will use quasar clustering observations as additional 
tests on our quasar models, as has been done in recent studies 
(e.g., Hopkins et al. 2008; Shankar et al. 2009a; Thacker et al. 
2009; BonoH et al. 2009; Croton 2009). 

The luminosity dependence of quasar bias can be modeled 
as follows. Assume at redshift z the probability distribution 
of host halo mass given bolometric luminosity L is p{M\L,z), 
then the halo averaged bias factor is: 

bdL,z) = J bM{M,z)p{M\L,z)dM , (3) 

where bmiM^z) is the linear bias factor for halos with mass 
M at redshift z. The derivation of the probability distribu- 



tion p{M\L,z) = dP{M\L,z)/dM is described in later sections 
[Eqn. (20)]. 

1 .3. Dark Matter Halo Mergers 

In the hierarchical universe small density perturbations 
grow under gravitational forces and collapse to form virial- 
ized halos (mostly consisting of dark matter). Smaller ha- 
los coalesce and merge into larger ones. Early work on 
the merger history of dark matter halos followed the ex- 
tended Press-Schechter (EPS) theory (e.g.. Press & Schechter 
1974; Bond et al. 1991; Lacey & Cole 1993), which is based 
on the spherical collapse model (Gunn & Gott 1972) with 
a constant barrier for the collapse threshold (the critical 
linear overdensity 6c « 1.68 is independent on mass, al- 
beit it depends slightly on cosmology in the ACDM uni- 
verse). Although the (unconditional) halo mass function 
n{M,z) predicted by the EPS theory agrees with numeri- 
cal A^-body simulations reasonably well, it is well known 
that it overpredicts (underpredicts) the halo abundance at 
the low (high) mass end (Sheth & Tormen 1999, and refer- 
ences therein). Motivated by the fact that halo collapses are 
generally triaxial, Sheth & Tormen (1999), based on earlier 
work by Bond & Myers (1996), proposed the ellipsoidal col- 
lapse model with a moving barrier where the collapse thresh- 
old also depends on mass. By imposing this mass depen- 
dence on collapse barrier, the abundance of small halos is 
suppressed and fitting formulae for the (unconditional) halo 
mass function are obtained, which agree with A^-body simu- 
lations much better than the spherical EPS predictions (e.g., 
Sheth & Tormen 1999; Sheth et al. 2001; Jenkins et al. 2001; 
Warren et al. 2006; Tinker et al. 2008). 

In addition to the unconditional halo mass function, the 
conditional mass function «(Mi,zi |Mo,zo) gives the mass 
spectrum of progenitor halos at an earlier redshift zi of a 
descendant halo Mq at redshift zo- This conditional mass 
function thus contains information about the merger history 
of individual halos and can be used to generate halo merger 
trees. A simple analytical form for «(Mi,zi |Mo,zo) can be 
obtained in the spherical EPS framework (e.g., Lacey & Cole 
1993); for the ellipsoidal collapse model, an exact but com- 
putationally consuming solution of the conditional mass func- 
tion can be obtained by solving the integral equation proposed 
by Zhang & Hui (2006). Recently, Zhang et al. (2008) de- 
rived analytical formulae for «(Mi,zi |Mo,zo) in the limit of 
small look-back times for the ellipsoidal collapse model. 

Alternatively, the halo merger rate can be retrieved directly 
from large volume, high resolution cosmological A^-body 
simulations, which bypasses the inconsistencies between the 
spherical EPS theory and numerical simulations. Using the 
product of the Millennium Simulation (Springel et al. 2005b), 
Fakhouri & Ma (2008) quantified the mean halo merger rate 
per halo for a wide range of descendant (at z = 0) halo mass 
10'^ <Mo< 1O'^M0, progenitor mass ratio lO"-' ^ C < 1 ^nd 
redshift < z < 6, and they provided an almost universal fit- 
ting function for the mean halo merger rate to an accuracy 
< 20% within the numerical simulation results: 



B{Mo,tz) 



Mo 



«(Mo,z) ' °'°^^^V1.2x10'2Mq 



0.083 



^-2.01 ^ 



exp 



0.098 



0.409 



V dz 



0.371 



(4) 



Here B(Mq,S^,z) is the instantaneous merger rate at redshift 
z, for halos with mass (Mo,Mo + iiMo) and with progeni- 
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tor mass ratio in the range (■^,■^ + '^0 [B(Mo,^,z) is in units 
of number of mergers x Mpc'^MQ^dz'^dt^], n(Mo,z) is the 
halo mass function (in units of Mpc"''Mq"'), ^ =M2/Mi < 1 
where Mi > M2 are the masses of the two progenitors, and 
Sc.z = ^c/D{z) is the Hnear density threshold for spherical col- 
lapse with D{z) the linear growth factor We adopt Eqn. (4) 
to estimate the halo merger rate in our modeling. The mass 
definition used here is friends-of-friends mass Mfof with a link 
length b = 0.2. For the range of mass ratios we are interested 
in (i.e., major mergers), the halo merger rate derived from the 
spherical EPS model can overpredict the merger rate by up to 
a factor of ~ 2 (Fakhouri & Ma 2008). 

1 .4. Galaxy ( subhalo ) Mergers 

Galaxies reside in the central region of dark matter halos 
where the potential well is deep and gas can cool to form stars. 
When a secondary halo merges with a host halo, it (along with 
its central galaxy) becomes a satellite within the virial radius 
of the host halo. It will take a dynamical friction time for 
the subhalo to sink to the center of the primary halo, where 
the two galaxies merge. The subsequent galaxy (subhalo) 
merger rate is therefore different from the halo merger rate 
discussed in §1.3, and a full treatment with all the dynamics 
(tidal stripping and gravitational shocking of subhalos) and 
baryonic physics is rather complicated. 

The simplest argument for the galaxy (subhalo) merger 
timescale is given by dynamical friction (Chandrasekhar 
1943; Binney & Tremaine 1987): 

/df ©orb A/host 
Tdf ~ — 7 — Tdyn , (J) 

In A Msat 

where Mhost and Msat are the masses for the host and satellite 
halos respectively, InA « ln(Mhost/A/sat) is the Coulomb log- 
arithm, 0oib is a function of the orbital energy and angular 
momentum of the satellite, f^f is an adjustable parameter and 
Tdyn = r/Ve{r) is the halo dynamical timescale, usually esti- 
mated at the halo virial radius ryir. This formula (5) is valid 
at the small satellite mass limit Mhost/A^sat ^ 1; although it 
is also used in cases of Mhost Msat in many semi-analytical 
models (SAMs) with the replacement of In A = (l/2)ln[l + 
(A^host/A^sat)^] (i-e., the original definition of the Coulomb log- 
arithm) or In A = ln(l +Mhost/Msat). However, in recent years 
deviations from the predictions by Eqn. (5) in numerical sim- 
ulations have been reported for both the Msat ^ Mhost and the 
Msat S Mhost regimes (e.g., Taffoni et al. 2003; Monaco et al. 
2007; Boylan-Kolchin et al. 2008; Jiang et al. 2008). In par- 
ticular even in the regime Msat /Mhost ^ 1 where the original 
Chandrasekhar formula is supposed to work, Eqn. (5) sub- 
stantially underestimates the merger timescale from simula- 
tions by a factor of a few (getting worse at lower mass ratios). 
This is because Eqn. (5) is derived by treating the satellite as 
a rigid body, while in reality the satellite halo loses mass via 
tidal stripping^^ and hence the duration of dynamical friction 
is greatly extended. 

Given the limitations of analytical treatments, we therefore 
retreat to numerical simulation results. Fitting formulae for 
the galaxy merger timescales within merged halos have been 
provided by several groups (Taffoni et al. 2003; Monaco et al. 
2007; Boylan-Kolchin et al. 2008; Jiang et al. 2008), and the 
merger rate of subhalos has also been directly measured from 
simulations (Wetzel et al. 2009; Stewart et al. 2008). 

^ The galaxy associated with the subhalo, on the other hand, does not suffer 
from mass stripping since it sits in the core region of the subhalo. 
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Fig. 1. — Merger timescales and (sub)halo merger rates. Upper, com- 
parison of the subhalo merger timescale from three dynamical friction pre- 
scriptions. Bottom: the redshift evolution of the halo merger rate (black solid 
lines), and the subhalo merger rate using the Boylan-Kolchin et al. (2008) 
(black dashed lines), and Jiang et al. (20()8) (red dashed lines) prescriptions 
for the merger timescale. for three (post)merger halo masses. 

The fitting formula for Tmerger in Jiang et al. (2008), who 
used hydro/A^-body simulations, has the following form: 

^_ 0.94r;°'^ + 0.6 1 _ 

where 77 is the circularity parameter [77 = (1 -e^)'/^, and e is 
the orbit eccentricity]. In deriving this equation Jiang et al. 
have removed the energy dependence of the orbit, i.e., they 
fix r^ = Tvir where r^ is the circular orbit that has the same 
energy as the satellite's orbit. In their simulation the ratio 
of r^/rvh ranges from ^ 0.6- 1.5 and has a median value ^ 
1. The distribution of circularity rj in their simulation (see 
their eqn. 7) has a median value 0.5, so in what follows 
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we take rj = 0.5. The redshift dependence of Tmergei is thus 
the same as that of the halo dynamical time Tdyn = rvii /Vf oc 

[Avi,-(z)]"'''2(l +zr^/^ (see Appendix A). 

Alternatively, the fitting formula of 
Boylan-Kolchin et al. (2008) is given by: 



' merger 



in 



,(e,z) = 0.216- 



exp(1.9?/) r^(E) _ 



"^dyn 



(7) 



Combining the halo merger rate (4) and galaxy (subhalo) 
merger timescale (6) or (7) we derive the instantaneous 
merger rate of galaxies within a merged halo of mass Mq for 
progenitor mass ratio ^ and at redshift z: 

figai(Mo,^,z) = B[Mo,^,Ze(z,0]$^ , (8) 

dz 

where z^iz, is a function of (z, and is determined by: 

^age(^) ^age(2e) ~ ^merger(C3^e) 7 (9) 

where fage(z) is the cosmic time at redshift z- Again, here 
Bgai(Mo,^,z) is the number of mergers per volume per halo 
mass per mass ratio per redshift. We find that c/ze/<^z is almost 
constant at all redshifts for ^ = 0.1 - 1 and monotonically de- 
creases as ^ increases. This constant can be approximated by 
its asymptotic value at large z when VL{z) — > 1 : 



dzst 
"5F 



— « <^ 1+0.22 



^ln(l + l/a 



In 2/3 



(10) 



for the fitting formula of Tmerger in Jiang et al. (2008), or 



dZe 

dz 



« <^ 1+0.09 



^'■^ln(l + l/0 



1 N 2/3 



(11) 



for the fitting formula of Tmeraer in Boylan-Kolchin et al. 
(2008). 

Unfortunately, current studies on the galaxy merger 
timescale have not converged yet, and different groups do re- 
port similar but quantitatively different results, at least partly 
caused by the different definitions and setups in their nu- 
merical simulations (e.g., hydro/A^-body versus pure A^-body 
simulations, halo finding algorithms and definition of merg- 
ers, etc). These fitting formulae are generally good within 
a factor of ~ 2 uncertainty. This uncertainty in the galaxy 
merger timescale within DM halos will lead to quite substan- 
tial differences in the galaxy merger rate at high redshift. To 
see this, we show in Fig. 1 the comparison of different fit- 
ting formulae for r,nerger and their consequences. In the up- 
per panel of Fig. 1 we show the ratio of Tmerger/Tdyn for the 
most likely orbit with circularity rj = 0.5 and rc = rvir using 
the fitting formula in Boylan-Kolchin et al. (2008) (solid line), 
Jiang et al. (2008) (dashed line), and a commonly-used SAM 
prescription: Tmerger = l-1677°'^'*[^ln(l + l/0]"'r(;/rvir (dotted 
line). For the major merger regime ^ > 0.3 the fitting for- 
mula in Boylan-Kolchin et al. (2008) agrees with the SAM 
prescription well, and they both approach the dynamical time 
at the high mass ratio end. However, the fitting formula in 
Jiang et al. (2008) yields a factor of ^ 2 longer than dynami- 
cal time at the high mass ratio end. 

In the lower panel of Fig. 1 we show the halo and galaxy 
major merger rates per unit time (integrated over ^ > 0.3), 
as function of redshift. The black solid lines show the 
halo major merger rate, the black dashed lines show the 
galaxy merger rate using the fitting formula of r^erger from 
Boylan-Kolchin et al. (2008), and the red dashed lines show 



the galaxy merger rate adopting the fitting formula from 
Jiang et al. (2008). In general the galaxy merger rates fall 
below the halo merger rate at high redshift and take over at 
lower redshift. At redshift z > 4, the galaxy merger rate us- 
ing the Jiang et al. (2008) formula is lower by almost an order 
of magnitude than that using the Boylan-Kolchin et al. (2008) 
formula (as well as the SAM prescription, since it agrees with 
eqn. 7 for ^ > 0.3), which makes it difficult to produce the 
quasar abundance at high redshift (see later sections). 

It is worth noting that although the merger rate peaks at 
some redshift, this trend is hierarchical such that it peaks at 
later time for more massive halos. This is somewhat in con- 
tradiction to the downsizing of the QSO luminosity function. 
This apparent discrepancy is likely caused by the fact that at 
low reshift, it becomes progressively more difficult for the 
black hole to accrete efficiently in massive halos because of 
the global deficit of cold gas due to previous mergers expe- 
rienced by massive halos and/or possible feedback quenches 
(e.g., Kauffmann & Haehnelt 2000; Croton et al. 2006). We 
will treat this fraction of QSO-triggering merger event as 
function of redshift and halo mass explicitly in our modeling. 

If we choose to normalize the total merger rate B (Bgai) by 
n{Mo,z), the abundance of descendant halos with mass Mq, 
at the same redshift z for halo mergers and galaxy mergers, 
then Bg^\/n falls below B/n at high redshift, then catches up 
and exceeds B/n a little, and evolves more or less parallel to 
B/n afterwards. Therefore the evolution in Bgai/n is shallower 
than that in B/n, consistent with the findings by Wetzel et al. 
(2009) once the different definitions of B/n and halo mass are 
taken into account. 

To summarize §1.3 and §1.4, we have compared the halo 
merger rate and galaxy merger rate as functions of redshift, 
(post)merger halo mass and mass ratio with different prescrip- 
tions for the dynamical friction timescale. It is, however, un- 
clear on which rate we should link to the QSO triggering rate. 
It is true that it will take a dynamical friction time for the two 
galaxies to merge after their host halo merged. But the black 
hole accretion might be triggered very early during their first 
orbital crossing, well before the galaxies merge. In what fol- 
lows, we adopt the halo merger rate to model the QSO trig- 
gering rate, and we discuss the consequences of adopting the 
alternative subhalo merger rates in §4.2. 

1 .5. BH Fueling During Mergers 

Modeling black hole fueling during mergers from first prin- 
ciples is not trivial: aside from the lack of physical inputs, 
most hydrodynamical simulations do not yet have the dynam- 
ical range to even resolve the outer Bondi radius. Although 
semi-analytical treatments of BH accretion are possible (e.g., 
Granato et al. 2004; Monaco et al. 2007), we do not consider 
such treatments in this paper because of the poorly understood 
accretion physics. 

Throughout this paper we adopt the following ansatv 

• QSO activity is triggered by a gas-rich major merger 
event (^^in < C < Ij- 

• The build-up of the cosmic SMBH population is mainly 
through gas accretion during the QSO phase. 

• The QSO light curve follows a universal form 
^'(ipeak, 0> where the peak luminosity Lpeak is correlated 
with the mass of the merged halo Mq. 
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TABLE 1 
Notation and Model Parameters 



Symbol 


Description 


ValueAJnits 


s 


Halo mass ratio 


0<5< 1 


Mo 


Postmerger halo mass 




M. 


Black hole mass 




SjBgal 


Halo(subhalo) merger rate 


Mpc-^Mg'rfr'dr' 


3>(L,z) 


QSO luminosity function 


Mpc-'dL-' 


Bl 


QSO triggeiing rate 


Mpc"VL"'f/z"' 


^eak 


QSO peak luminosity 






QSO triggeiing redshift 




/qSo(J,Mo) 


QSO triggering fraction 




The Lpcak-'Wo relation 


7 


Slope 


5/3 


Cfe = 0) 


Nonnalization at z = 


log(6 X 10''5)-127 




Scatter 


0.28 dex 


Ca) = C(z = 0) + A log(l+z) 


Evolution in normalization 


/3, =0.27=1/3 


The Hght curve model 



Radiative efficiency 0.1 

1 n38 ,^ 

fSalpctcr e-foldmg time 



/ Eddington luminosity per Mq 1.26 X lO"** ergs 'M^-' 



/(l-e)Ao 

/ Fraction of seed BH mass to peak BH mass 10 ^ 

?peak Time to reach the peak luminosity (- In /)?salpeter 

Afl Eddington ratio before fp^ak 3 

a Power-law slope of the decaying phase 2.5 

The QSO-triggering rate 

Minimum mass ratio 0.25 

Mn,in Exponential lower mass cut 3 X lO" /i^'Mq 

Mmax(z) = Mq„ench(l +^)'' ■■ Exponential Upper mass cut 10'^(1 A^'Mg 



Note. — Parameter values are for the fiducial model. 

Once we have specified the QSO light curve model, and es- 
timated the QSO triggering rate from the major merger rate, 
we can convolve them to derive the QSO luminosity func- 
tion. Within this evolutionary QSO framework, we can derive 
the distributions of host halo masses and BH masses for any 
given instantaneous luminosity and redshift, allowing us to 
compare with observations of quasar clustering and Edding- 
ton ratio distributions. The details of this framework are elab- 
orated in §2. 

We note that, any correlations established with our model 
are only for hosts where a gas-rich major merger is triggered 
and self-regulated BH growth occurs. Host halos which expe- 
rience many minor mergers or dry mergers may deviate from 
these relations. We will come back to this point in the discus- 
sion section. 

2. MODEL FORMALISM 

In this section we describe our model setup in detail. Some 
notations and model parameters are summarized in Table 1 
for clarity. 

2.1. Determining the Luminosity Function 

Major mergers are generally defined as events with mass 
ratio ^tnin = 0.3 < ^ < 1, but our framework can be easily 
generalized to other values of ^niin- Because it takes more 
than just a major merger event to trigger QSO activity, we 
shall model the QSO triggering rate as a redshift and mass- 
dependent fraction of halo merger rate. The QSO luminosity 
function at a given redshift z can be obtained by combining 



the triggering rate and the evolution of QSO luminosity : 

^(L,z)dL= I BL[Lp,,i,(L,t,-t,,),z']dLp,,kdz' , (12) 

where ^{L,z)dL is the QSO number density in the luminos- 
ity range (L,L + dL) at redshift z; Bi{L^eskiZ') is the QSO- 
triggering rate (merger number density per redshift per peak 
luminosity) at redshift z' with peak bolometric luminosity 
Lpesk{L,t^-t,i), which is determined by the specific form of 
the light curve; and t^i are the cosmic time at z and z' ■ In 
integrating this equation we impose an upper limit of redshift 
Zmax = 20, but we found that the integral converges well be- 
low this redshift since the merger rate decays rapidly at high 
redshift. Eqn. (12) implies that the distribution of triggering 
redshift z' given L at redshift z is: 

dP{L^^,z'\L,z) , , dLpesk 

=p(Lpe„k,Z |L,Z) CxBL(ipeak,Z ) , 

(13) 

where again Lpeak is tied to L via the light curve model. Obvi- 
ously only quasars with Lpeak > L can contribute to $(L,z). 
The QSO-triggering rate BiiL-p^^^^z) is related to the halo 

Note that we have implicitly assumed that a second QSO-triggering 
merger event does not occur during the major accretion phase of the QSO 
- a reasonable assumption since statistically speaking bright QSOs are short- 
lived (fgso <S Ih)- Observationally binary /multiple quasars within a single 
halo with comparable luminosities are rare occurrences (/binary ^ 0.1%, e.g., 
Hennawi et al. 2006, 2009; Myers et al. 2008), further supporting this as- 
sumption. 
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merger rate B(Mo,^,z) by: 

fiL(iiDeak,z)t/i'peak =/qSo(z,Mo) / B{Mo,S,,z)d^dMQ 

^min 

where we parameterize the fraction /qso as 

Mniin(z) Mo 



(14) 



/Qso(z,Mo) = :^(z)exp 



Mo M„,x(z), 



(15) 



where jr(2) describes how /qso(z,Mo) decreases towards 
lower redshift, i.e., the overall reduction due to the consump- 
tion of cold gas. We also introduce exponential cutoffs of 
/qso at both high and low mass such that at each redshift, 
halos with too small Mo cannot trigger QSO activity; on the 
other hand, overly massive halos cannot cool gas efficiently 
and BH growth halts, and the gas-rich fraction may become 
progressively smaller at higher halo masses (especially at low 
redshift). We discuss choices for these parameters later in §3. 

To proceed further we must specify the relation between 
ipeak and Mo, and choose a light curve model. We assume 
that the Lpeak-Mo correlation is a power-law with log-normal 
scatter, as motivated by some analytical arguments and hy- 
drodynamical simulations (e.g., Wyithe & Loeb 2002, 2003; 
Springel et al. 2005b; Hopkins et al. 2005; Lidz et al. 2006): 



fifP(LpeaklMo) 
t/logLpeak 

= (2wi)-i/2 



exp 

where the mean relation is: 

^peak 



= /?(ipeak|Mo) 

[log (C+7logMo)]- 



2al 



log 



ergs 



= C+7log 



Mo 
/i-'M, 







,(16) 



(17) 



with normalization C and power-law slope 7. The log scatter 
around this mean relation is denoted as ct/,. When the scatter 
between Lpeak and Mo is incorporated, Eqn. (14) becomes: 



■Sl (ipeak, 2) 

/qso 



B(Mo,tz)di 



Mo 

^peak 



/:'(Lpeak|Mo)t/logMo , 



(18) 



from which we obtain the probability distribution of post- 
merger halo mass Mo at fixed peak luminosity Lpeak and red- 
shift z: 



dP(Mo\Lp^^i,,z) 
d log Mo 



/7(Mo|Lpeak,z) 

Mo 



/qso / B{Mo,^,z)d^ 



-'peak 



p(Lpeak|Mo) . (19) 



We derive the probability distribution of host halo mass at 
given instantaneous luminosity and redshift, p{Mo\L,z), as: 



dP(Mo\L,z) 



= p{Mo\L,z) 



(i log Mo 

£/P(Mo|Lpeak,z') dP(L^,,i„z'\L,z) 



c/ log Mo 



dz' 



dz! , (20) 



which can be convolved with the halo linear bias bM{Mo,z) to 
derive the QSO linear bias biiL, z) at instantaneous luminosity 
L and redshift z- 



Since we have assumed that a second QSO-triggering 
merger event has not occurred, at redshift z the postmerger 
halo Mo should maintain most of its identity as it formed 
sometime earlier, i.e., we neglect mass added to the halo via 
minor mergers or accretion of diffuse dark matter between the 
halo formation time and the time when the QSO is observed 
at redshift z- 

Given the halo mass distribution (20), we can further derive 
the "active" halo mass function hosting a QSO with L > Lmin: 



d^ 



Mo 



dlogMo 



HL,z)dL 



dP{Mo\L,z) 
dlogMo 



(21) 



which can be used to compute the halo duty cycles. 

2.2. The Light Curve Model 

For the light curve model there are several common 
choices in merger-based cosmological QSO models (e.g., 
Kauffmann & Haehnelt 2000; Wyithe & Loeb 2002, 2003): 
a) a light bulb model in which the QSO shines at a constant 
value Lpeak for a fixed time fQso; b) an exponential decay 
model L = Lpeak exp(-r/f qso); c) a power-law decay model 
L = Lpeak(l +f/fQso)~" with a > 0. Once a light curve model 
is chosen, the total accreted mass during the QSO phase is 
simply^ : 



_ 1-e 

M, j-elic ~ 9~ 



L{Lp^iik,t)dt . 



(22) 



The evolution of the luminosity Eddington ratio A = L/L^^^ 
is (neglecting the seed BH mass): 



A(f) = 



L(LpeakjO _ 

lM,{t) ~ 1(1 -e) 



L(Lpealci t) 



(23) 



where / = 1.26 x 10^'^ergs~'Mg is the Eddington luminosity 
per Mq . 

Unfortunately none of the three LC models is self- 
consistent without having A 1 at the very beginning of 
BH growth, simply because they neglect the rising part of 
the light curve. To remedy this, we must model the evo- 
lution of luminosity and BH growth self-consistently (e.g.. 
Small & Blandford 1992; Yu & Lu 2004, 2008). We consider 
a general form of light curve where the BH first grows ex- 
ponentially at constant luminosity Eddington ratio Ao (e.g., 
Salpeter 1964) to Lpeak at f = fpeak, and then the luminosity 
decays monotonically as a power-law (e.g., Yu & Lu 2008): 



Lpeak exp 



L(Lpeak,0= < 



— J—^{t-tpeak) 



-^peak 



^peak / 

where fpeak is determined by: 
Lpeak = /AoM,,o exp 



,0<f <fpeak 

(24) 

^ ^ ^peak ; 



/(l-e)Ao^ 



"^peak 



(25) 



^ Throughout the paper we assume constant radiative efficiency e. At 
the very late stage of evolution or under certain circumstance (i.e., hot gas 
accretion within a massive halo), a SMBH may accrete via radiatively- 
inefficient accretion flows (RIAFs) with very low e and mass accretion rate 
(e.g., Narayan & Yi 1995). We neglect this possible accretion state since the 
mass accreted during this state is most likely negligible (e.g. Hopkins et al. 
2006). 



8 



SHEN 



where M, o is the seed BH mass at the triggering time f = 
0. In eqn. (24), a determines how rapidly the LC decays. 
Larger a values lead to more rapid decay. Thus this model 
accommodates a broad range of possible light curves. With 
this Ught curve model (24) we have: 



(L, 



M.(Lpeak,0^ 



peak 

"7a7 



exp 



^ (r-fpeak; 



^ ^peak ^ ( 1 ^)^peak^peak ^ 



, < f < fpeak 



(26) 



1 



l-a 



^peak / 



■1 



t > t, 



peak 



where we require a > 1 so that a BH cannot grow infinite 
mass. The evolution of the luminosity Eddington ratio is 
therefore: 



< f < f , 



peak 



Mt)=< 



(f /fpeak) " 



1 . 1(1 -e)t, 



peak 



Ao ec^(l-a) 





1-Q -1 


\(-) 


-1 


. \ fpeak / 





peak 



(27) 



where during the decaying part of the LC the Eddington ratio 
monotonically decreases to zero. We assume that the seed BH 
mass is a fraction / of the peak mass M,(fpeak) = i'peak/CAo), 
and we have 



fpeak ~ " 



l{l-e)Xo 



ln/ = (-ln/)fsalpeter 



(28) 



In what follows we set / = 10 . Thus the seed BH is negligi- 
ble compared to the total mass accreted, and fpeak ~ 6.9fsaipetei- 
where fsaipeter is the (e-folding) Salpeter timescale (Salpeter 
1964). Although Ao = 1 is the formal definition of Eddington- 
limited accretion, it assumes the electron scattering cross sec- 
tion, and in practice super-Eddington accretion cannot be 
ruled out (e.g., Begelman 2002). So we consider Ao within 
the range log Ao G [-1, 1]. 
The relic BH mass, given eqn. (26), is simply 



. - _ ^peak (1 ^)^peakfpeak _ ^eak 

"*»,relic — -7T I Z 77 5 — "77 

/Ao (a-l)ec- /Ao 



1- 



In/ 

a-l 



(29) 



If instead we choose an exponential model for the decay- 
ing part of the LC (e.g., Kauffmann & Haehnelt 2000), L{t > 
fpeak) = ipeak exp [-(f- fpeak) /fQSo], then the rclic mass be- 
comes: 

Lpenk _|_ (1 ~ £)^eakfQSO 



relic = 



/Ao 



(30) 



Therefore LC models with large values of a mimic the expo- 
nentially decaying model, and we do not consider the expo- 
nential LC model further. 

If the decaying parameter a is independent of mass, then 
Eqns. (26) and (29) imply M, i-eiic oc Lpeak- In other words, 
the scalings between Lpeak and Mo, and between M. ^iic and 
Mo should be the same for QSO-triggering merger remnants 
(e.g., early type galaxies). Some theoretical models and simu- 
lations predict Lpeak k M^^^ (e.g.. King 2003; Di Matteo et al. 
2005; Springel et al. 2005a), where momentum is conserved 
in self-regulated feedback; while others invoking energy 

conservation predict Lpeak ocMq'^ (e.g.. Silk & Rees 1998; 
Wyithe & Loeb 2003). Two recent determinations of the 
local black hole mass-halo mass relation gave M. ieiic oc 




0.0 0.5 



1.0 1.5 2.0 

Time (Gyr) 



2.5 3.0 



Fig. 2. — Example light curves for e = 0. 1 and Ao = 3.0. for a = 1.5, 2.5, 3.5. 
Larger values of a lead to more rapid decay. The time that the QSO spends 
above half of its peak luminosity is short, < 100 Myr for large values of a. 

Ml'^^ (Ferrarese 2002) and M.jeiic oc M^ -^ (Baes et al. 2003). 
Within uncertainties M, relic oc Lpeak seems to be a reasonable 
scaling, therefore we do not consider further mass dependence 
of a. Given the definition of fpeak (Eqn. 28), our light curve 
model (24) is then universal in the sense that it scales with 
Lpeak in a self-similar fashion. 

As an example. Fig. 2 shows three light curves with a = 
1.5,2.5,3.5 for e = 0.1 and Ao = 3.0, in which case fpeak ^ 
115 Myr. The time that a QSO spends above 50% of its 
peak luminosity is typically < 100 Myr for sharp decaying 
curves, but it can be substantially longer for extended decay- 
ing curves. 

2.3. BH Mass and Eddington Ratio Distributions 

Now we have specified the light curve model and the peak 
luminosity-halo mass relation, and tied the luminosity func- 
tion $(L,z) to the QSO-triggering merger rate. We continue 
to derive the instantaneous BH mass and Eddington ratio dis- 
tributions at fixed instantaneous luminosity L and redshift 
z, which can be compared with observationally determined 
distributions (e.g., Babic et al. 2007; Kollmeier et al. 2006; 
Shen et al. 2008a; Gavignaud et al. 2008). 

Suppose that at redshift z we observe quasars with instan- 
taneous luminosity L. These quasars consist of objects trig- 
gered at different earlier redshift z' and are at different stages 
of their evolution when witnessed at z (e.g., Eqn. 12). There 
is a characteristic earlier redshift Zc, determined by 

fage(z)-fage(Zc) = fpeak , (31) 

where fage is the cosmic time. Quasars triggered between 
[Zc,z] are all in the rising part of the LC; while quasars trig- 
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-1 



logL (ergs ) 



Fig. 3. — Predicted bolometric luminosity functions z) = /dlogL in our fiducial model (solid lines). Overplotted are the complied bolometric LF data 
from Hopkins et al. (2007) for optical (black circles), soft X-ray (green squares) and hard X-ray (red triangles) samples. The dashed lines are the predicted LF 
for a model where the normalization in the mean Lpeak~Mo relation is reduced by an additional amount at i > 3.5, tuned to fit the LF at z > 4.5 (see §4.2 for 
details). The blue dotted lines are the predictions based on an alternative prescription for the redshift evolution in the normalization of the Lpeak —Mq relation as 
discussed in §4.4. 



gered earlier than Zc are all in the decaying part of the LC. In 
order to contribute to <i?(L,z), a quasar should have peak lu- 
minosity Lpeak = if it is triggered right at Zc, and it should 
have Lpeak > L otherwise. Therefore the triggering redshift 
distribution (iP(Lpeak,2'|i',2)/i^2' (Eqn. 13) peaks around^ z' « 
Zc since BLC^peakiZ') decreases rapidly when Lpeak increases. 
Moreover, at the observing redshift z, all the contributing 
quasars triggered between [Zc,z] will have Eddington ratio 
Xi -'^- = \q and BH mass M,x.;'^; = L/{1\q); while all the 
contributing quasars triggered before Zc will have Eddington 
ratio Al,z'^z < Ao and BH mass M,x.z'^z > L/{lXo). The 
probability distribution of instantaneous BH mass M, at the 

* Except for low luminosities when the lower halo mass cut Mmin in the 
QSO-triggering rate starts to kick in, the distribution dP(L^.^^,z'\L,z)/dz' 
instead has a dip around Zc- 



observing redshift z and luminosity L is therefore: 
dP dz! 
dP(M,\L,z) _ d^d\ogM, 



dlogM, 



1- 



' dP 
. d^ 



dz 



M. > — , (32) 



where dP(z'\L,z)/dz' is given by Eqn. (13) and dz! /dM, is 
determined by the decaying half of the LC (Eqn. 27), i.e., the 
triggering redshift z'{> Zc) is a monotonically increasing func- 
tion of M, (note that here M, refers to the instantaneous BH 
mass observed at z). The denominator in the above equation is 
to normalize the distribution - we literally augment the distri- 
bution with the pileup of objects with M, = L/(/Ao) triggered 
within [ZctZ]- The fraction of such objects is about 20% for 
L > 10"*^ ergs"' and becomes negligible at lower luminosities. 
This procedure removes the spike (S function) at the constant 
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BH mass M, = L/(/Ao) in the distribution, which is an artifact 
of our moder. Eqn. (32) gives the equivalent distribution of 
Eddington ratios at z and at instantaneous luminosity L. 

Given the luminosity function (12) and the distribution of 
M, (32) we can determine the active BH mass function (i.e., 
the mass function of active QSOs above some luminosity cut 

^min) 



M. 



dlogM, 



<^>{L,z)dL 



dP{M,\L,z) 
dlogM, 



(33) 



3. THE REFERENCE MODEL 



We are now ready to convolve the LC model with the QSO- 
triggering rate Eqn. (18) to predict the LF ^(L,z) using Eqn. 
(12). Before we continue, let us review the model parameters 
and available observational/theoretical constraints. 

• QSO triggering rate ACmin, M„,i„{z), M^^Azl ^iz)]: We 
choose our fiducial value of ^min ~ 0.3 following the 
traditional definition of major mergers (1 : 3 or 1 : 4). 
For the minimum halo mass Mmin below which efficient 
accretion onto the BH is hampered we simply choose 
Mniin(z) = 3 X lO"/r'M0. We model the maximum 
halo mass cut above which the gas-rich merger fraction 
drops as a function of redshift; 



Mmax(z)=Mquench(l+Z> 



(34) 



with /3 > 0. Therefore at lower redshift halos have a 
smaller upper threshold. Since at z < 0.5 rich group 
to cluster size halos no longer host luminous QSOs we 
tentatively set Mquench = 1 x 10'^/!~'Mq. The parame- 
ters Mmin and M^ax control the the shape of the LF at 
both the faint and bright luminosity ends. The global 
gas-rich merger fraction, !F{z), is more challenging to 
determine. Direct observations of the cold gas frac- 
tion as function of redshift and stellar mass (and there- 
fore halo mass) are still limited by large uncertainties. 
Moreover, how gas-rich a merger needs to be in order 
to trigger efficient BH accretion is not well constrained. 
For these reasons, we simply model T{z) as a two-piece 
function such that J^{z) = 1 when z>2, the peak of the 
bright quasar population; and J^{z) linearly decreases to 
J-{) at z = 0. We note again that /qso(2,-^o) should be 
regarded as the fraction of major mergers that trigger 
QSO activity. 

• The Lpeak-Mo relation (C, 7, (Jl): some merger event 
simulations (e.g., Springel et al. 2005a; Hopkins et al. 
2005; Lidz et al. 2006) reveal a correlation between the 
peak luminosity and the mass of the (postmerger) halo: 
(Vak) ~ 3 X 1045(Mo/10i2riMo)4/3ergs-' at z = 2, 
with a lognormal scatter ct/, = 0.35 dex. Other analyt- 
ical arguments predict 7 = 5/3, where energy is con- 
served during BH feedback (e.g.. Silk & Rees 1998; 
Wyithe & Loeb 2003). We choose 7 values between 
[4/3,5/3] since the above simulations are also consis- 
tent with the 7 = 5/3 slope. Furthermore, we allow 
the normalization parameter C to evolve with redshift, 
C(z) = C{z = 0) + log[(l +z)A] with /3i > 0. Thus for 

^ We have tested with the alternative BH mass distribution at instanta- 
neous luminosity L and redshift ; where the contribution of objects with 
M, = L/(1\q) is described by a 5 function with normalization J~ ^dz! , and 
we find that other derived distributions based on this BH mass distribution 
are almost identical to those using Eqn. (32). 



the same peak luminosity, the host halos become less 
massive at higher redshift. This decrease of the charac- 
teristic halo mass with redshift is also modeled by sev- 
eral other authors, although we do not restrict to /3i = 1 
(e.g., Lapi et al. 2006), or more complicated prescrip- 
tions (e.g., Wyithe & Loeb 2003; Croton 2009), since 
the form of the L-Mq (or M, -Mq) relation is also dif- 
ferent in various prescriptions. A smaller characteris- 
tic halo mass is easier to account for the QSO abun- 
dance at high redshift since there are more mergers of 
smaller halos, but it also reduces the clustering strength. 
The intrinsic scatter of the Lpeak-Mo relation, ct/,, will 
have effects on both the luminosity function and clus- 
tering. Larger values of <tl lead to higher QSO counts, 
but will also dilute the clustering strength due to the up- 
scattering of lower mass halos (e.g.. White et al. 2008). 
The Lpeak-Mo relation establishes a baseline for the 
mapping from halos to SMBHs. Although in our fidu- 
cial model we adopt the above rather simple scaling 
[oc (1 +z)^'] for the redshift evolution in C(z), we will 
discuss alternative parameterizations for C(z) in §4.2 
and §4.4. 

• The light curve model ( Xq, e, fpeak> 01} : our chosen fidu- 
cial value for fpeak is fpeak = 6.9fSalpeter (/ = 10"^); but 

our results are insensitive to the exact value of /. The 
radiative efficiency is e « 0.1 from the Soltan argument 
(Eqn. 1). We choose Ao between [0.1, 10] and a > 1.1. 

Normally to find the best model parameters one needs to 
perform a minimization between model predictions and 
observations. We do not perform such exercise here because 
it is difficult to assign relative weights to different sets of ob- 
servations (i.e., LF, quasar clustering, Eddington ratio distri- 
butions, etc), and we don't know well enough the systematics 
involved. Instead, we experiment with varying the model pa- 
rameters within reasonable ranges, to achieve a global "good" 
(as judged by eye) fit to the overall observations. Our fidu- 
cial model has the following parameter values: ^min = 0.25, 
jTo = LO, /3 = 1.5, 7 = 5/3, C(z = 0) = log(6 x lO'*^)- I27, 
(3i =0.27, (Tl = 0.28, Ao = 3.0 and a = 2.5, which are all within 
reasonable ranges. In particular we found that the parame- 
terization of JF(z) is unnecessary because the effect of cold 
gas consumption is somewhat described by Mjnax(z) already. 
Below we describe in detail the predictions of this reference 
model and comparison with observations, and we defer the 
model variants and caveats to §4. 

Fig. 3 shows the model LF d'^/dlogL = Lln(10)<I>(L,z) in 
solid lines, given by Eqn. (12), where we overplot the com- 
piled bolometric luminosity function data in Hopkins et al. 
(2007). In computing Eqn. (12) we integrate up to z = 20 
but the integral converges well before that. The model 
under-predicts the counts at the faint luminosity end (L < 
10"*^ ergs"') for z < 0.5, leaving room for low luminosity 
AGNs triggered by mechanisms other than a major merger 
event (i.e., by secular processes). At redshift z > 4.5, the 
model also underestimates the QSO abundance. This could be 
alleviated if the characteristic halo mass shifts to even lower 
values, i.e., even larger values of the parameter C at z > 4.5 
(see discussions in §4.2 and §4.4). Alternatively, the high-z 
SMBH population may be different in the sense that it is not 
tied to ^ > 0.3 major mergers events directly. Also, the halo 
merger rate (Eqn. 4) has been extrapolated to such high red- 
shifts for the massive halos considered here. Despite these 
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Fig. 4. — Distributions of triggering redsliift (upper panels) and lialo mass (lower panels) in our fiducial model, for thi'ee instantaneous luminosities and two 
values of observing redshifts. The distributions of triggering redshift generally peak around Zc (marked by arrows) given by Eqn. (31). For lower luminosity, the 
distribution of the triggering redshift is more extended. As a consequence, there is a significant population of faded low-luminosity QSOs which have massive 
hosts. 



0.5 



10'*'*ergs-i at z > 2.5 



facts, the model correctly reproduces the LP from z 
to z « 4.5. The turnover below L 

is caused by the lower mass cutoff Mmin = 3 X 10" /i-'Mq. At 
lower redshift, this turnover flattens due to the gradual pile up 
of evolved high-peak luminosity quasars well after fpeak- 

We show some of the predicted distributions in Fig. 4 for 
this reference model. The upper panels show the distributions 
of the triggering redshift z' for instantaneous luminosities 
L= 10'*^10'^^10'^^ergs-' atz= 1 (upper left) andz = 2(upper 
right). As expected, the distribution of the triggering redshift 
peaks around the characteristic redshift Zc given by Eqn. (31). 
The bottom panels show the distributions of host halo mass 
Mo for the three instantaneous luminosities, at the two red- 
shifts respectively. For bright quasars (L > lO'^^ergs"'), the 
distributions of Mq are roughly log-normal, with the width 



and mean determined mainly through the Lpeak-A/o relation 
- but also slightly modified through the convolutions with 
other distributions (see Eqns. 13, 19 and 20). For faint QSOs 
(L < lO'^^ergs"'), however, the halo mass distribution has a 
broad high-mass tail contributed by evolved high-peak lumi- 
nosity quasars triggered earlier This contamination of mas- 
sive halos hosting low luminosity QSOs will increase the clus- 
tering bias of the faint quasar population. 

Fig. 5 shows our model predictions for the red- 
shift and luminosity dependence of quasar bias, com- 
pared with observations. The model-predicted red- 
shift evolution of quasar bias is in good agreement 
with observations (e.g., Porciani et al. 2004; Croom et al. 
2005; Porciani & Norberg 2006; Myers et al. 2006, 2007a,b; 
Shen et al. 2007; da Angela et al. 2008; Padmanabhan et al. 
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Fig. 5. — Predictions for the lineal' bias of quasar clustering in our fiducial 
model. Upper, bias evolution with redshift. The data points are from the 
measurements in Shen et al. (2009) and three lines show the predicted bias 
for quasars with instantaneous luminosity L= Iff*-^ , 10"**, 10"*' ergs"' . The 
median luminosity of quasars used in the clustering analysis for the six red- 
shift bins are: (log(L/ergs"')) = 45.6,46.3,46.6,46.9,46.9,47.1. Bottom: 
predicted luminosity dependence of quasar bias. The filled circles are from 
the measurements in Shen et al. (2009) for the 10% most luminous and the 
remainder of a sample of quasars spanning 0.4 < z < 2.5. The open cir- 
cles are the measurement of quasar clustering for 2.9 < z < 3.5 in Shen et al. 
(2009). The open square is from Francke et al. (2008) using cross-correlation 
of X-ray selected AGN with high redshift galaxies, and the open triangles ai'e 
from Adelberger & Steidel (2005) using cross-correlation of optical AGNs 
with galaxies. Note that for the quasar clustering data we plot biases derived 
from both with and without including negative correlation function data in 
the fitting (see table 1 of Shen et al. 2009). 

2008; Shen et al. 2008b, 2009; Ross et al. 2009). But it has 
some difficulties in accounting for the large bias at z « 4, re- 
sulting from the need to reproduce the quasar abundance at 
such high redshift (i.e., the characteristic halo mass shifts to 
lower values). We discuss possible solutions to this problem 
in §4.2. Our model also predicts the luminosity dependence 
of quasar clustering. At low luminosities, quasar bias de- 
pends weakly on luminosity, which is consistent with obser- 
vations (e.g., Porciani & Norberg 2006; Myers et al. 2007a; 
da Angela et al. 2008; Shen et al. 2009). This simply re- 
flects the fact that quasars are not light bulbs, i.e., there 
is non-negligible scatter around the mean Lpeak-M) rela- 
tion, and some faint quasars live in massive halos because 
of the evolving light curve (e.g., Adelberger & Steidel 2005; 
Hopkins et al. 2005; Lidz et al. 2006). On the other hand, 
the model predicts that the bias increases rapidly towards 
high luminosities and/or at high redshift. This is consistent 
with the findings by Shen et al. (2009) that the most luminous 
quasars cluster more strongly than intermediate luminosity 
quasars at z < 3. Our model-predicted luminosity dependent 
quasar clustering broadly agrees with the measurements in 



Shen et al. (2009) for 0.4 < z < 2.5, but we caution that their 
measurements are done for samples with a broad redshift and 
luminosity range in order to build up statistics (e.g., see their 
fig. 2), hence a direct comparison is somewhat difficult. At 
redshift ~ 3, Adelberger & Steidel (2005) and Francke et al. 
(2008) measured the clustering of low luminosity AGN us- 
ing cross-correlation with galaxies. Their data are shown 
in Fig. 5 in open square and triangles for the measurements 
in Francke et al. (2008) and Adelberger & Steidel (2005) re- 
spectively, which are broadly consistent with our predic- 
tions. However the clustering of optically bright quasars 
is apparently at odds with the low bias value derived from 
the AGN-galaxy cross-correlation in Adelberger & Steidel 
(2005), likely caused by the fact that the latter AGN sam- 
ple spans a wide redshift range 1-6 < z < 3.7 and luminosity 
range, and that the uncertainty in their bias determination is 
large. 

Fig. 6 shows our model predictions for the Eddington ratio 
distributions at various redshifts and instantaneous luminosi- 
ties. Aside from the cutoff at Ao = 3 (our model setup), these 
Eddington ratio distributions are approximately log-normal, 
and broaden towards fainter luminosities. This again reflects 
the nature of the light curve - at lower luminosities, there 
are more objects at their late evolutionary stages and shin- 
ing at lower Eddington ratios. These predictions are in good 
agreement with the observed Eddington ratio distributions 
for bright quasars (L > lO'^^ergs"') where the distribution 
is narrow and peaks at high mean values ((log A) G [-1,0]) 
(e.g., Woo&Urry 2002; Kollmeier et al. 2006; Shenetal. 
2008a), as well as for faint AGNs (L < lO'^^ergs"') where 
the distribution is broader and peaks at lower mean values 
((log A) e [-3,-1]) (e.g., Babicetal. 2007; Gavignaud et al. 
2008). However the predicted mean values and widths are 
not necessarily exactly the same as those determined from 
observations, since uncertainties in the BH mass estimators 
used in these observations may introduce additional scatter 
and biases in the observed Eddington ratio distributions (e.g., 
Shen et al. 2008a). The Eddington ratio distributions of faint 
AGNs (L < 10"*"* ergs"') are particularly broad at low redshift 
since more high-Lpeak objects have had enough time to evolve. 
On the other hand, at high redshift z > 2, the minimum allow- 
able Eddington ratio is set by the cosmic age at that redshift. 



e.g., A„ 



10 at z = 2 using Eqn. (27). This explains the 



narrowing of the lower end of the Eddington ratio distribution 
forL= !()"*'* ergs"', seen in the right panel of Fig. 6. 

Since we have reproduced the observed luminosity function 
and predicted the Eddington ratio distributions as function of 
luminosity, we can use Eqn. (33) to derive the active BHMF 
in QSOs above some minimum luminosity. In Fig. 7 we show 
the total BHMF (setting L^in = in Eqn. 33) assembled at 
several redshifts, where the gray shaded region indicates the 
estimates of the local dormant BHMF based on various galaxy 
bulge-BH scaling relations (Shankaret al. 2009b). Our model 
prediction for the local BHMF is incomplete at the low-mass 
end, mainly because our model does not include low lumi- 
nosity AGN activity (presumably linked to M, < lO^M© BH 
growth) possibly triggered by secular processes (see further 
discussion in §4.2), and also because we set the minimal 
halo mass that can trigger QSO activity during major merg- 
ers Mmin = 3 X 10" /i"'Mq in order to reproduce the faint-end 
LF at high redshift. At the high mass end, the predicted slope 
in the total BHMF broadly agrees with (although is shallower 
than) that for the local dormant BHMF. The majority of the 
present-day > 10^-^ Mq BHs were already in place by z = 1, 
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Fig. 6. — Predictions for the Eddington ratio distributions in our fiducial model for three different redshifts. The distribution is narrower and has higher peak 
values towards higher luminosities. These distributions are in good agreement with the observed distributions for both bright quasars and low luminosity AGNs 
(e.g., Kollmeier et al. 2006; Babic et al. 2007; Shen et al. 2008a; Gavignaud et al. 2008). 



but less than 50% of them were assembled by z = 2. This 
is somewhat in disagreement with McLure & Dunlop (2004), 
who claimed that the majority of > 10^-^ SMBHs are al- 
ready in place at z 2 based on virial BH mass estimates of 
optically selected bright quasars. We suspect this discrepancy 
is caused by: 1) the fact that virial BH mass estimates tend 
to systematically overestimate the true BH masses due to a 
Malmquist-type bias (see discussions in Shen et al. 2008a); 
2) the high-mass end slope in the local dormant BHMF is 
steeper than our model predictions. On the other hand, our 
model predictions are in good agreement with the predictions 
by Shankar et al. (2009b) for the high-mass end. It is possi- 
ble to make our model prediction agree with the local BHMF 
at the high-mass end by imposing some cutoff in the light 
curve [Eqn. (24)] such that BHs cannot grow too massive; but 
a more accurate observational determination of the high-mass 
end of the local BHMF is needed to resolve these issues. 

Fig. 8 shows the halo duty cycles, defined as the ratio of 
the number density of active halos hosting QSOs brighter 
than Lmin to that of all halos, as function of halo mass. We 
have used the Sheth et al. (2001) halo mass function and Eqn. 
(21) for the active halo mass function. For quasars (L > 
10^^ ergs"') and for typical halo mass Mo - 2 X lO'^/r'Mo, 
the duty cycle is ~ 0.15,0.1,0.03,0.01 at z = 3,2, 1,0.5. 

We can also compare the model predicted BH mass distri- 
butions for quasars within certain luminosity ranges with the 
virial BH mass estimates. For this purpose we use the virial 
BH mass estimates from Shen et al. (2008a) for the SDSS 
DR5 quasar catalog (Schneider et al. 2007). These optical 
quasars are flux-limited to / = 19.1 at z < 3 and we have used 
eqn. (1) in Shen et al. (2009) to convert /-band magnitude to 
bolometric luminosity. For quasars/AGNs where the SMBHs 
are still actively accreting, what we observe is the instanta- 
neous BH mass rather than the relic BH mass. In Fig. 9 
the solid lines show the model predictions for the BH mass 
distributions of SDSS quasars at z = 1 and z = 2, weighted 
by the LF; the dashed lines show the distributions based on 
virial BH masses. It is remarkable that not only the distribu- 
tions of our model predictions are broader, but also the peaks 
are shifted to lower masses (~ 0.6 dex) compared with those 
from virial mass estimates - as already discussed extensively 
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Fig. 7. — BHMFs assembled at various redshifts in our fiducial model. 
The gray shaded region shows the estimates for the local BHMF from 
Shankar et al. (2009b). The red hne shows the prediction for the local BHMF, 
which is incomplete at M, < Iff ^ Mq by a factor of a few because we did 
not include contributions from AGNs triggered by secular processes or minor 
mergers (as reflected in the failure to reproduce the LF at the low luminos- 
ity end L < 10'*' ergs"' at z < 0.5, see Fig. 3). The yellow and green lines 
show the predicted local BHMF at M. > 10** /i^'Mq after correcting for BH 
coalescence; these corrections are likely upper limits (see §4.3 for details). 

in Shen et al. (2008a, i.e., the Maknquist-type bias in virial 
mass estimates). 

4. DISCUSSION 

4.1. Comparison with Previous Work 

The merger basis of our framework, as advocated by many 
authors (e.g., Kauffmann & Haehnelt 2000; Wyithe & Loeb 
2002, 2003; Volonteri et al. 2003), provides a physical origin 
for the QSO population, and distinguishes the current study 
from other works which focus on BH growth using the QSO 
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Fig. 8. — Halo duty cycles, i.e., fraction of hales hosting QSOs brighter than Ln,i„, 
Sheth et al. (2001) and the active halo mass function from Eqn. (21). 



as function of halo mass, computed using the halo mass function from 



LF as an input (e.g., Yu & Tremaine 2002; Yu & Lu 2004, 
2008; Marconi et al. 2004; Merloni 2004; Shankar et al. 2004, 
2009b). Our simple framework, although semi-analytical 
in nature, accommodates a wide range of updated and new 
observations of QSO statistics; these new observations in- 
clude quasar clustering and Eddington ratio distributions. 
Most of the early quasar models (e.g., Haiman & Loeb 1998; 
Kauffmann & Haehnelt 2000; Wyithe & Loeb 2002, 2003; 
Volonteri et al. 2003) focused mainly on the luminosity func- 
tion (partly because other observations were not available at 
that time), and most of them assumed simplified light curve 
models which were unable to reproduce the observed Edding- 
ton ratio distributions. 

Hopkins et al. (2008) presented a merger-based quasar 
model that utilizes a variety of quasar observations for com- 
parison, including the latest clustering measurements. Our 
model framework is different from theirs in two major as- 



pects: 1) they estimated the major- merger rate from the com- 
bination of empirical halo occupation models of galaxies 
and merging timescale analysis, while we used directly the 
halo merger rate from simulations - both approaches have 
their own advantages and disadvantages; 2) the light curve in 
their model is extracted from their merger-event simulations 
(Hopkins et al. 2005), while we have adopted a parametriza- 
tion of the light curve which is fit by observations. Their light 
curve model is clearly more physically-motivated than ours, 
yet it needs to be confirmed in future simulations with higher 
resolution and better understandings of BH accretion physics. 
On the other hand, the virtue of our framework is that it allows 
fast and easy estimations of model predictions and parameter 
adjustments to fit updated observations. 

4.2. Caveats in the Reference Model 
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Fig. 9. — Distributions of BH masses in quasars at two redshifts. Tlie solid 
lines are the predictions of our reference model for i < 19.1 (i.e., the flux 
limit in the main SDSS quasar catalog) quasars. The dashed lines show the 
distributions of virial BH masses from Shen et al. (2008a). 

As already mentioned in §3, our fiducial model is not an 
actual fit to the overall observations (LP, clustering, and 
Eddington ratio distributions) because of the ambiguity of as- 
signing relative weights to individual observational data sets. 
Instead, we have experimented with varying the model pa- 
rameters within reasonable ranges, to achieve a global "good" 
(as judged by eye) fit to observations. If consider only the 
LP as observational constraints, there are model degenera- 
cies between the luminosity decaying rate a and the nor- 
malization of the mean Lpeak-A^o relation C, i.e., if QSO 
luminosity decays more slowly, the typical host halos need 
to shift to more massive (and less abundant) halos in order 
not to overpredict the LP; likewise, if the scatter around the 
mean Lpeak-A^o relation increases, we can reduce C to 
match the LP. However, these degeneracies are broken once 
the clustering observations are taken into account. A large 
scatter <tl or slow varying light curve (small values of a) 
cannot fit the large bias at high redshift and the luminos- 
ity dependence of clustering. Our model is also more com- 
plicated than previous models (e.g., Haiman & Loeb 1998; 
Kauffmann & Haehnelt 2000; Wyithe & Loeb 2002, 2003; 
Volonteri et al. 2003; Shankar et al. 2009b) in the sense that 
we have more parameters, which is required for a flexible 
enough framework to accommodate a variety of observations. 

Por our fiducial model we have used the halo merger rate 
as the proxy for the QSO-triggering rate. Alternatively, if 
we use the subhalo merger rate (the delayed version of halo 
merger rate; §1.4), it makes little difference below z = 3, but 
it further underestimates the LP at z > 3, as expected from 
Pig. 1. The bolometric LP data we used here has uncertain- 
ties both from measurements and bolometric corrections at 
the ~ 20-30% level (cf., Hopkins et al. 2007; Shankar et al. 
2009b), not enough to reconcile the discrepancy at z > 4.5. 
There are several ways to modify our model to match ob- 
servations of LP at z > 4.5: a) decrease the threshold ^,„jn 
above which mergers can trigger QSO activity; b) modify the 
i'peak-^o relation such that at fixed peak luminosity, QSOs 



shift to even smaller halos at z > 4.5, which can be achieved 
by including some higher order terms in the redshift evolu- 
tion of C; and c) increase the scatter so that more abun- 
dant low mass halos can contribute to the LP by up-scattering. 
As an example of such a model, we modify the redshift evo- 
lution of the mean Lpeak-A/o relation at high redshift such 
that C C + log[(i^f)^/2] ^ > 3 5 -^j^j^ ^^^^ additional 
term of evolution in the mean Lpeak - relation, the typical 
host halo shifts to lower masses (Mq 4 x 1O"/!~'M0 for 
(logL/ergs"') = 46 at z = 6, compared to Mq ~ 1O'^/i"'M0 in 
our fiducial model), and the resulting LP (dashed lines in Fig. 
3) fit the observations well. This modification has little effects 
on the LP at z < 3.5, as well as other predicted QSO proper- 
ties. An alternative parametrization of the redshift evolution 
in C is further discussed in §4.4. 

However, there are no independent constraints on the mass 
of halos hosting quasars at z > 4.5 (such as those inferred 
from quasar clustering), and the merger scenario of quasar 
activity may be different at such high redshift. Hence we do 
not attempt to fully resolve this issue in this paper 

There is also slight tension between quasar clustering and 
LP at z ~ 4 in our model, as already noted by several stud- 
ies (White et al. 2008; Wyithe & Loeb 2009; Shankar et al. 
2009a). On the one hand, we need smaller halo masses to 
account for quasar abundance; we also need larger halos to 
account for the strong clustering on the other To fully re- 
solve this issue we need better understanding of the halo bias, 
and reliable fitting formulae for it, derived from simulations 
for the relevant mass and redshift ranges'*, as well as better 
measurements of quasar clustering at high redshift with fu- 
ture larger samples. Nevertheless, our fiducial model is still 
consistent with both LP and clustering observations within the 
errors. 

Our model also underpredicts the LP at the low lumi- 
nosity end [L < 10"*^ ergs"') at z < 0.5. This is somewhat 
expected, since our model does not include contributions 
from AGNs triggered by mechanisms other than a major 
merger. The fuel budget needed to feed a low luminosity 
AGN (L < lO'^^ergs"') is much less stringent than that for 
bright quasars. Therefore secular processes (i.e., gas inflows 
driven by bars, tidal encounters, stochastic accretion, as well 
as minor mergers), while not as violent and efficient as major 
mergers, provide viable means to fuel AGNs at low activity 
levels. At z < 0.5, there are evidence that some low luminos- 
ity AGN [L 10'*-'~'*"'ergs"', powered by intermediate-mass 
BHs M, ~ 10^ Mq) hosts have no classical bulges'^, which 
indicates that secular processes are responsible for trigger- 
ing these low-luminosity AGN activity and BH growth (e.g., 
Greene et al. 2008, and references therein). There are also 

* For instance, it has been suggested that in addition to mass, halo cluster- 
ing also depends on concentration, assembly history and recent merger activ- 
ity (the so-called assembly bias, e.g., Gao et al. 2005; Wechsler et al. 2006; 
Wetzel et al. 2007). If recently merged halos (as quasar hosts) have larger 
bias than average for the same halo mass, then it may reconcile the slight ten- 
sion between reproducing both the LF and quasar clustering at high redshift 
(e.g., Wyithe & Loeb 2009), although current estimate of this enhancement 
is only on the level of ~ 10% (e.g., Wetzel et al. 2007). 

' Classical bulges (including ellipticals and bulges in early-type disk galax- 
ies) are presumably formed via mergers and they follow the fundamental 
plane of elhptical galaxies (e.g.. Bender et al. 1992). On the other hand, sec- 
ular processes can build up the so-called pseudobulges, which have distinct 
structural properties from classical bulges (Kormendy & Kennicutt 2004). 
While pseudobulges can host central BHs, there are some evidence that the 
bulge-BH mass relations for pseudobulge systems are offset from that for 
classical bulge systems (e.g., Hu 2008; Greene et al. 2008), indicating that 
secular processes are less efficient in building BHs. 



16 



SHEN 



observational implications that the bulk of BH growth has 
shifted from the most massive BHs (M, > 10^ Mq) at high 
redshift (z < 2) to low mass BHs (M. < 10^ Mq) locally (e.g., 
McLure & Dunlop 2004; Heckman et al. 2004), along with 
the cosmic downsizing in luminosity function evolution (e.g., 
Steffenetal. 2003; Uedaetal. 2003; Hasinger et al. 2005; 
Hopkins et al. 2007; Bongiorno et al. 2007). The typical 
transition from merger-driven BH growth to secularly-driven 
BH growth likely occurs around BH mass ~ lO'' Mq (e.g., 
Hopkins & Hernquist 2009), corresponding to L ~ 10 ergs"' 
at the Eddington limit. Since the specific merger rate in- 
creases towards higher redshift, while the rate of secular pro- 
cesses is almost constant with time and these processes are 
relatively slow and inefficient for BH growth, it is conceiv- 
able that secular processes will only become important at late 
times (z < 0.5 for instance) in building up the low mass end 
of the SMBH population. The implementation of secularly- 
driven AGN activity in our model will be presented in future 
work. 

Finally, we mention that our model predicts a turnover in 
the LF at L < 10"*^ ergs"' at z > 3. This simply reflects the 
lower mass cutoff at M^jn = 3 x 1O"/2"'M0 in our model. 
Future deeper surveys for low luminosity AGNs at z > 3 are 
necessary to probe this luminosity regime, and to impose con- 
straints on how efficiently SMBHs can form in low mass ha- 
los. 

4.3. The Effects ofBH Coalescence on the BHMF 

The BHMFs discussed in §3 and Fig. 7 are the mass func- 
tions from accretion only, i.e., we have neglected the effects of 
BH coalescence. As discussed in §1.1, BH coalescence will 
redistribute the BH mass function but does not change the to- 
tal BH mass density (neglecting mass loss via gravitational 
radiation) since essentially all the mass ended up in BHs were 
accreted. There are two routes in which the coalescence with 
pre-existing BHs might become important within our model 
framework'". First, both galaxies in the merging pair of ha- 
los are elliptical, i.e., they already experienced a QSO phase 
in the past and formed massive nuclear BHs. In this case the 
current merger event will be a dry merger and form a SMBH 
binary without triggering a new QSO. Second, during the ma- 
jor merger, one galaxy is elliptical and the other one is spiral, 
in which case a QSO will be triggered and the mass of the old 
BH of the previous elliptical will add to the new BH system. 

A complete exploration of these routes and their conse- 
quences (including the effects of BH ejection, mass loss 
through gravitational radiation, etc.) can be better achieved 
with Monte Carlo realizations of halo merger trees and our 
model prescriptions for QSO triggering and BH growth, 
which we plan to investigate in a future paper Here we can 
approximately estimate the maximal impact of BH coales- 
cence on the BHMFs in the two cases, assuming that BH coa- 
lescence always occurs within a Hubble time and all the mass 
in the pre-existing BHs is added to the final BH, i.e., no BH 
ejection or mass loss via gravitational radiation. We focus on 
the high-mass end (M, > IQ^Mq) of the local BHMF since 
our model prediction is incomplete at the lower-mass end. 

Although the dry merger case is not subject to the ma- 
jor merger condition for QSO-triggering, we still restrict to 

In contrary to the two cases discussed below, we assume that a major 
merger between two spirals will lead to negligible BH mass contribution from 
the pre-existing BHs, since in our model setting, spiral galaxy has not yet 
experienced a major merger and hence significant BH growth. 



C > Cmin here because in minor mergers: 1) it will take too 
long for the two pre-existing BHs to become a close binary, 
and 2) the mass increment due to BH coalescence is insignif- 
icant. Observational determination of the major dry merger 
rate is difficult, and the current best estimate is; on aver- 
age, present-day spheroidal galaxies with My < -20.5 (cor- 
responding to M, > 10** M0) have undergone 0.5-2 major 
dyr mergers since z ~ 0.7 (Bell et al. 2006). If we assume all 
the M, > 10^ M0 BHs undergo one 1 : 1 dry merger after the 
QSO phase, the local BHMF will redistribute as the yellow 
line in Fig. 7. In this case the abundance of the most massive 
(M, > a few X lO^M©) BHs is enhanced by up to a factor of 
~2 atM. = 10'°/i"'Mo. 

In the half-dry major merger case, the fraction of the pre- 
existing BH mass to the final BH mass is (1 + ^')~^'''' ~ 0.07- 
0.7, where 1/4 < ^' < 4 is the major merger mass ratio (in our 
reference model) between the two halos. This is because the 
final BH mass after a QSO phase scales as the 5/3 power to 
the halo mass in our model. Averaging over possible values 
of the fractional increment due to the pre-existing BH is 
~ 30%. If all QSO-triggering mergers are this kind of half- 
dry event, the predicted z = BHMF (red line in Fig. 7) will 
redistribute from lower mass to higher mass (the green line). 
The enhancement at the high-mass end of the local BHMF is 
comparable to the dry major merger case. In practice QSO- 
triggering mergers can occur between two spirals, hence the 
actual correction due to pre-existing BHs in half-dry mergers 
should be smaller 

Combining these two cases we conclude that the impact 
of BH coalescence on the BHMF is probably insignificant 
compared with other uncertainties and systematics in cur- 
rent observations and our model framework. Similar con- 
clusions were also achieved in several independent work 
(Volonteri et al. 2003; Yu & Lu 2008; Shankar et al. 2009b) 
albeit with difference in details. 

4.4. Implications for BH Scaling Relations 

In our formalism there is a generic scaling relation between 
the relic BH mass and halo mass, which has the same slope 
and scatter as the Lpeak-A^o relation (from Eqns. 17 and 29): 

The local M, - Mq relation for dormant BHs reported in 
Ferrarese (2002) and Baes et al. (2003) has a slope in the 
range 1.3-1.8, and a normalization lower by a factor of 
~ 3 - 5 than our predictions. This is probably due to the fact 
that we neglected continued growth of halos by minor merg- 
ers and diffuse matter accretion since the major merger event, 
if most of the local massive BHs were assembled at z > 1 (as 
our model predicts). The fact that we did not impose a cutoff 
in the LC (24) may also lead to overly massive BHs. We note 
that although the normalization (and perhaps slope as well) 
of the local M, -Mq relation may depend on galaxy morpho- 
logical type (e.g., Zasov et al. 2004; Courteau et al. 2007; Ho 
2007), early-type galaxies (SO and ellipticals) appear to oc- 
cupy the upper envelope in the M,-Mq relation. Our model 
scaling relations are only valid for early type galaxies which 
are presumably merger remnants. 

However, it should be pointed out that the BH-halo scaling 
relation relies on the assumed Lpeak-A^o relation. The simple 
prescription for its evolution in our reference model already 
shows some difficulties in reproducing the QSO LF at z > 4.5 



SMBHS IN THE HIERARCHICAL UNIVERSE 



17 



(Fig. 3; see §4.2). Hence our fiducial model is not very ap- 
propriate for predicting the redshift evolution of the BH-halo 
scaling relation. To make this point more clear, let us consider 
a more physically-motivated prescription for the Lpeak-A/o re- 
lation in which L oc V^^^, where Vvii is the halo virial velocity 
(Eqn. A5), and the normalization of the Lpeak-A^o relation 
evolves as (Wyithe & Loeb 2002, 2003): 



C(z) = C(z = 0) + ^log(l+z) + ^log 



Avir(z) 



Avu(0) 



(36) 



e.g., the evolution in C{z) is more rapid than our fiducial set- 
ting. With this new implementation for the L^^h^-Mq rela- 
tion, we found that a good global fit can be achieved with 
the following parameter adjustments (other parameters are the 
same as in our reference model): Mquench = 3 x 10'^/i~'Mq, 
C(z = 0) = log(0.8 X 10"*^)- 127 and crz. = 0.4(1 +z)-i/l Since 
the redshift evolution in the normalization C{z) is now faster, 
the starting value C(0) is reduced; consequently the expo- 
nential upper cut of halo mass, Mquench, increases in order to 
account for QSO counts at low redshift. The scatter in the 
i'peak - Mq relation needs a redshift evolution to achieve ad- 
equate fits for both clustering and LP over a wide redshift 
range. The predicted LP is shown as dotted lines in Pig. 3, 
which does a much better job at z > 4.5 than our fiducial 
model. Other predicted properties are slightly degraded (but 
still are reasonably good fits to observations) than those pred- 
icated by our fiducial model. 

This new L^^^-Mq relation, together with the approxima- 
tion that the halo virial velocity Vvii ~ Vc, the galaxy circu- 
lar velocity, and the assumption that the local Vc-cr relation 
(Perrarese 2002) does not evolve, result in a constant M, - a 
relation (Wyithe & Loeb 2003). Neglecting the scatter, the lo- 
cal mean Vc -cr relation in Perrarese (2002) is: 



log 



( 

\kms" 



= 0.84 log 



kms 



+ 0.55 



(37) 



Thus the new Lpeak-A/o relation predicts a constant M,-a 
relation (e.g., Eqns. 17, 26, 29, 36, and A6): 



M, 



10** 



:0.9~5.0x 



200 km s- 



4.2 



(38) 



with a slope and normalization (bounded by M, peak and 
M, relic) consistent with local estimates (Gebhardt et al. 2000; 
Perrarese & Merritt 2000; Tremaine et al. 2002; Lauer et al. 
2007). This consistency, however, is built on a couple of 
assumptions and approximations, and neglecting successive 
evolution after the self-regulation of BH and bulge growth. 
Given all these complications, it is beyond the scope of the 
current study to fully settle this issue. We simply remind the 
reader that a non-evolving M, - a relation is generally allowed 
within our framework. 

5. CONCLUSIONS 

We have developed a general cosmological framework for 
the growth and cosmic evolution of SMBHs in the hierarchi- 
cal merging scenario. Assuming that QSO activity is trig- 
gered by major mergers of host halos, and that the resulting 
light curve follows a universal form with its peak luminos- 
ity correlated with the (post)merger halo mass, we model the 



QSO LP and SMBH growth self-consistently across cosmic 
time. We tested our model against a variety of observations of 
SMBH statistics: the QSO luminosity function, quasar clus- 
tering, quasar/ AGN BH mass and Eddington ratio distribu- 
tions. A global good fit is achieved with reasonable parame- 
ters. We summarize our model specifics as follows: 

• The QSO-triggering rate is determined by the ^ > 0.25 
halo merger rate with exponential cutoffs at both the 
low and high halo mass ends M^in = 3 x 1O"/!~'M0, 
M^,Az)=W\l+zfl-h-'MQ. 

• The universal light curve follows an initial expo- 
nential Salpeter growth with constant Eddington ra- 
tio Aq = 3 for a few e-folding times to reach the 
peak luminosity Lpeak, which is correlated with the 
(post)merger halo mass as (Lpeak) = 6 x 10^*^(1 + 
z)'/^(Mo/lO'^/z"'M0)^/^ergs"\ with a log-normal 
scatter ctl = 0.28 dex. It then decays as a power-law 
with slope a = 2.5. 

Our simple model successfully reproduces the LP, quasar 
clustering, and Eddington ratio distributions of quasars and 
AGNs at 0.5 < z < 4.5, supporting the hypothesis that QSO 
activity is linked to major merger events within this redshift 
range. However, there are still many unsettled issues. Be- 
low we outline several possible improvements of our simple 
model, which will be addressed in future work. 

Our model under-predicts the LP at the faint luminosity end 
L < lO'^^ergs"' at z < 0.5, which is linked to the growth of 
the less massive < 10^ Mq SMBHs. This is indicative of a 
population of low luminosity AGNs triggered by mechanisms 
other than major mergers at the low redshift universe - either 
by minor mergers or secular processes. We need to incorpo- 
rate this ingredient in our SMBH model, in order to match the 
local BHMP at the low mass end (M, < 10^ M©). 

In our modeling we have neglected the possibilities of a 
closely-following second major merger event and the trig- 
gering of two simultaneous QSOs during a single major 
merger event. Therefore our model does not include more 
than one QSOs within a single halo. We will use Monte- 
Carlo realizations of halo merger trees to assess the proba- 
bility of such rare occurrences and see if they can account for 
the small (< 0.1%) binary /multiple quasar fraction observed 
(Hennawi et al. 2006, 2009; Myers et al. 2008). 

Our model can be improved to include the radio loudness 
of QSOs as well. If radio loudness requires both a massive 
host halo (to provide the hot IGM) and a massive SMBH (to 
launch the kinetic jet), we can statistically populate radio-loud 
QSOs in halos within our model framework. It can be tested 
against the clustering of radio-loud quasars (e.g., Shen et al. 
2009) and the radio-loud fraction as function of luminosity 
and redshift (e.g., Jiang et al. 2007, and references therein). 

We thank the anonymous referee, Prancesco Shankar, 
Michael Strauss, Scott Tremaine, and Martin White for 
constructive comments that have greatly improved the 
manuscript. We are grateful to Prancesco Shankar for point- 
ing out an error in the merger rate equation (4) in an earlier 
version of the manuscript. This work was supported by NSP 
grant AST-0707266. 
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APPENDIX 
DARK MATTER HALOS 

All dark matter halos are assumed to have a spherical NFW profile (Navarro et al. 1997): 



(r/r,)(l + r/r,)2 



where rj and ps are the characteristic scale and the density at this scale. 
The virial mass and virial radius are related by (Bullock et al. 2001): 



47r 



-AvirP„r; 



.3 

vir ' 



(Al) 



(A2) 



where Avir(z) ~ {Idiir^ + di2x - 39x^) / n{z) with x = ftiz)- 1 (Bryan & Norman 1998) is the spherical overdensity relative to the 
background matter density p,, and ri(z) = [1 +r2A(l +z)~^/Jlo]~'- Note there is the slight difference in defining the virial radius 
in Bullock et al. (2001) and Navarro et al. (1997) where the latter uses raoo (the radius corresponding to a spherical overdensity 
200 times the critical density) to define the virial radius. Both definitions of virial radius are frequently used in studies on the 
galaxy merger time scale within merged dark matter halos: Jiang et al. (2008), Stewart et al. (2008) used the former definition, 
while Boylan-Kolchin et al. (2008), Wetzel et al. (2009) used the latter. 
The enclosed mass within an NFW profile truncated at r„ is 



M(r„) : 



drAiir^ pnv^{r) = ^Tipsr] 



ln(l+c)- 



1+c 



where c = ro/r,. Therefore we have 



Mvii- = M(rvir) = 47rp,rj 



ln(l+Cvir)- 



l+c„ 



(A3) 



(A4) 



where Cvk = r^ijr^ is the usual definition of the concentration parameter. The mean relation between Myi, and Cvii is given in 
Bullock et al. (2001). 

The virial velocity (usually defined as the circular velocity at the virial radius) Vvh, and the maximum circular velocity at 
'"max ~ 2.16r5 are (Bullock et al. 2001): 



vir 



2 _T/2 



0.216c„ 



'•vir ' V4 ln(l + Cvir)-Cvir/(l+Cvii) 

This imphes that the relation between virial mass Mvii and virial velocity Kii is: 



Mvir(z) : 



Ait 



-Avir(z)p( 



-1/2 



(1 +zr^'^G-^I^Vl, = 1.37 X IO'^I^q'/VMq 



200kms- 



(1+z)" 



-3/2 



Avir(z) 
Avir(O) 



-1/2 



where po = 2.78 x lO^fio/i^A^oMpc ^ is the z = mean matter density. 
The DM halo dynamical time r^yn is usually defined as ryir/Vvir" 



Tdyn = = 1.4 X lO'^yr x [VL^h'^Al +z)1 

'vir 



31-1/2 



(A5) 



(A6) 



(A7) 



For simplicity we have neglected the difference between Mvii and the friends-of-friends mass Mfof (with a link length b = 0.2) 
throughout the paper. But we give an approximate conversion formula below for completeness: 

Mm ^ ln(l-l-Cfof)-Cfof/(l-l-Cfof) 

Mvir ln(l-|-Cvir)-Cvir/(l+Cvir) ' 

where cm = /"fof/r? can be solved via the following equation: 



(A8) 



Cfof(l+Cfof) 



9 ln(l-|-Cvir)-Cvir/(l+Cvir) ' 

where we have used the fact that the density at r^f is pNFw('"fof) ~ 3p„/(27r/7^). 



(A9) 
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